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^— ^ Abstract 

Fifty years ago, Fitzugh introduced a phase portrait that became famous for a twofold rea- 
^ son: it captured in a physiological way the qualitative behavior of Hodgkin-Huxley model and 

U it revealed the power of simple dynamical models to unfold complex firing patterns. To date, 

in spite of the enormous progresses in qualitative and quantitative neural modeling, this phase 
portrait has remained the core picture of neuronal excitability. Yet, a major difference between 
^ the neurophysiology of 1961 and of 2011 is the recognition of the prominent role of calcium 

channels in firing mechanisms. We show that including this extra current in Hodgkin-Huxley 
I— I dynamics leads to a revision of Fitzugh-Nagumo phase portrait that affects in a fundamental 

way the reduced modeling of neural excitability. The revisited model considerably enlarges 
the modeling power of the original one. In particular, it captures essential electrophysiological 
• signatures that otherwise require non-physiological alteration or considerable complexication 

of the classical model. As a basic illustration, the new model is shown to highlight a core 
dynamical mechanism by which the calcium conductance controls the two distinct firing modes 
of thalamocortical neurons. 

^ 1 Introduction 

> 

Rooted in the seminal work of Hodgkin and Huxley [H], conductance-based models have become 
a central paradigm to describe the electrical behavior of neurons. These models combine a num- 
ber of advantages, including physiological interpretability (parameters have a precise experimental 
meaning) and modularity (additional ionic currents and/or spatial effects are easily incorporated 
using the interconnection laws of electrical circuits [TUlf5]). Not surprisingly, the gain in quantitative 
description is achieved at the expense of mathematical complexity. The dimension of detailed quan- 
titative models makes them mathematically intractable for analysis and numerically intractable for 
the simulation of large neuronal populations. For this reason, reduced modeling of conductance- 
based models has proven an indispensable complement to quantitative modeling. In particular, the 
FitzHugh-Nagumo model [5] , a two-dimensional reduction of Hodgkin-Huxley model, has played an 
essential role over the years to explain the mechanisms of neuronal excitability (see e.g. [521 [5] for an 
excellent introduction and further references). More recently, Izhikevich has enriched the value of 
reduced-models by providing the Fitzugh-Nagumo model with a reset mechanism |13j that captures 
the fast (almost discontinuous) behavior of spiking neurons. Such models are used to reproduce 
the qualitative |S1I1] and quantitative [HlUT] behavior of a large family of neuron types. Notably, 
their computational economy makes them good candidates for large-scale simulations of neuronal 
populations 

The Hodgkin-Huxley model and all reduced models derived from it H] focus on sodium and 
potassium currents, as the main players in the generation of action potentials: sodium is a fast 
depolarizing current, while potassium is slower and hyperpolarizing. Initally motivated by reduced 
modeling of dopaminergic neurons in which calcium currents are essential to the firing mechanisms 
[7] , the present paper mimicks the classical reduction of the Hodgkin-Huxley model augmented with 
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an additional calcium current. The calcium current is a distinct player because it is depolarizing, 
as the sodium current, but acts on the slower timescale of the potassium current. 

To our surprise, the inclusion of calcium currents in the HH model before its planar reduction 
leads to a novel phase portrait that has been disregarded to date. Mimicking earlier classical work, 
we perform a normal form reduction of the global HH reduced planar model. The mathematical 
normal form reduction is fundamentally different in the classical and new phase portrait because it 
involves a different bifurcation. The classical fold bifurcation is replaced by a transcritical bifurca- 
tion. 

The results of these mathematical analysis lead to a novel simple model that further enriches 
the modeling power of the popular hybrid model of Izhikevich. A single parameter in the new 
model controls the neuron calcium conductance. In low calcium conductance mode, the model 
captures the standard behavior of earlier models. But in high calcium conductance mode, the 
same model captures the electrophysiological signature of neurons with a high density of calcium 
channels, in agreement with many experimental observations. For this reason, the novel reduced 
model is particularly relevant to understand the firing mechanisms of neurons that switch from a 
low calcium-conductance mode to a high calcium-conductance mode. Because thalamocortical (TC) 
neurons provide a prominent example of such neurons, they are chosen as the main illustration of 
the present paper, the benefits of of which should extend to a much broader class of neurons. 

2 Planar reduction of Hodgkin-Huxley model revisited in 
the light of calcium channels 

Calcium channels participate in the spiking pattern by providing, together with sodium channels, 
a source of depolarizing currents. In contrast to sodium channels whose gating kinetics are fast, 
calcium channels activate on a slower time-scale, similar to that of potassium channels [TT]. As a 
consequence, they oppose the hyperpolarizing effect of potassium current activation, resulting in 
bidirectional modulation capabilities of the post-spike refractory period. We model this important 
physiological feature by considering the HH model [12] with an additional voltage-gated (non- 
inactivating) calcium current Ica smd a DC-current Ipump that accounts for hyperpolarizing calcium 
pump currents. The resulting model is similar to HH model when the conductance of the calcium 
currents is low, but becomes strikingly different when the calcium conductance is high. Note that 
the inactivation of calcium channels is not included in the HH dynamics because it is known to take 
place in a much slower time scale [25] . The inactivation will typically be accounted for by a slower 
adaptation of the calcium conductance, see Section [5] 

Figure la illustrates the spiking behavior induced by the action of an external square current 
lapp in the two different modes. As compared to the original HH model (Fig. la left), the presence 
of the calcium current is characterized by a triple electrophysiological signature (see Fig. la right): 

• spike latency: the spike train (burst) is delayed with respect to the onset of the stimulation 

• plateau oscillations: the spike train oscillations occur at a more depolarized voltage than the 
hyperpolarized state 

• after-depolarization potential (ADP): the burst terminates with a small depolarization 

This electrophysiological signature is typical of neurons with sufficiently strong calcium currents. See 
for instance: spike latency [20, J^, plateau oscillations [2 , ADPs 1, 5 . However, the mechanisms 
by which these behaviors occur have never been analyzed using reduced planar models to date. 

Following the standard reduction of HH model [5], we concentrate on the voltage variable V 
(that accounts for the membrane potential) and on a recovery variable n (that accounts for the 
overall gating of the ion channels) as key variables governing excitability (see methods). The 
phase-portrait of the reduced HH model is shown in Figure lb (left). This phase portrait and 
the associated reduced dynamics are well studied in the literature (see [S] for the FitzHugh paper, 
and [HI H] for a recent discussion and more references). We recall them for comparison purposes 
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Figure 1: Step responses of the HH model without (left) and with a calcium current (right). 

(a) Time-evolution of the applied excitatory current (bottom) and of the corresponding membrane potential 
(top) in HH model (the reduced model leads to almost the same behavior (Supplementary Fig. SI)), (b) 
Phase portraits of the reduced Hodgkin-Huxley model in resting (top) and spiking states (bottom). The 
V- and n-nuUclines are drawn as a full and a dashed line, respectively. Trajectories are drawn as solid 
oriented red lines. Black circles denote stable fixed points, white circles unstable fixed points, and cross 
saddle points. The presence of calcium channels strongly affects the phase-portrait and the corresponding 
electrophysiological time-response of the neuron to excitatory inputs. 
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only. The resting state is a stable focus, which lies near the minimum of the familiar N-shaped 
y-nuUcline. When the stimulation is turned on (spiking mode), this fixed point loses stability in 
a subcritical Andronov-Hopf bifurcation (see the numerical bifurcation analysis of Supplementary 



Section S.2), and the trajectory rapidly converges to the periodic spiking limit cycle attractor. As 



the stimulation is turned off (resting mode), the resting state recovers its global attractivity via a 
saddle-nodc of limit cycles (the unstable one being born in the subcritical Hopf bifurcation), and 
the burst terminates with small subthreshold oscillations (cf. Fig. la left). 

In the presence of the calcium current, the phase-portrait changes drastically, as shown in Figure 
lb (right). In the resting mode, the hyperpolarized state is a stable node lying on the far left of 
the phase-plane. The y-nuUcline exhibits a "hourglass" shape. Its left branch is attractive and 
guides the relaxation toward the resting state after a single spike generation. The sign of V changes 
from positive to negative approximately at the funnel of the hourglass, corresponding to the ADP 
apex. The right branch is repulsive and its two intersections with the n-nullcline are a saddle and 
an unstable focus. 

When the stimulation is turned on, the y-nullcline breaks down in an upper and a lower branch. 
The upper branch exhibits the familiar N-shape and contains an unstable focus surrounded by a 
stable limit cycle, very much as in the reduced Hodgkin-Huxley model. In contrast, the lower branch 
of the y-nullcline, which is not physiological without the calcium currents, comes into play. While 
converging toward the spiking limit cycle attractor from the initial resting state, the trajectory must 
travel between the two nuUclines where the vector field has smaller amplitude. As a consequence, 
the first spike is fired with a latency with respect to the onset of the stimulation, as observed in 
Figures la (right) in the presence of the calcium current (see also Supplementary Fig. ^I]). In 
addition, a comparison of the relative position of the resting state and the spiking limit cycle in 
Figure lb (right) explains the presence of plateau oscillations. As the stimulation is turned off the 



spiking limit cycle disappears in a saddle- homoclinic bifurcation (see Supplementary Sections S.2 
and S.3), and the resting state recovers its attractivity. 

The presence of the lower branch of the y-nullcline has a physiological interpretation. In the 
reduced HH model, the gating variable n accounts for the activation of potassium channels and the 
inactivation of sodium channels. Their synergy results in a total ionic current that is monotonically 
increasing with n for a fixed value of V (Supplementary Fig. ^left). In this situation, at most one 
value of n solves the equation V = and there can be only one branch for the voltage nullcline. 
In contrast, when calcium channels are present, the reduced gating variable must capture two an- 
tagonistic effects. As a result, the total ionic current is decreasing for low n (the gating variable is 
excitatory), and increasing for large n (the gating variable recovers its inhibitory nature) (Supple- 
mentary Fig. fright). In this situation, two distinct values of n solve the equation V = 0, which 
explains physiologically the second branch of the ^-nullcline. To summarize, the lower branch of 
the voltage nullcline accounts for the existence of an excitatory effect of n, which is brought by 
calcium channel activation. 



3 The central ruler of excitability is a transcritical bifurca- 
tion, not a fold one 

The power of mathematical analysis of the reduced planar model ^ is fully revealed by introducing 
two further simplifications. 

• Time-scale separation: we exploit that the voltage dynamics are much faster than the recovery 
dynamics by assuming a small ratio n ~ 0{e)V (the approximation holds away from the 
voltage nullcline) and by studying the singular limit e = 0. 

• Transcritical singularity: by comparing the shape of the voltage nullcline in Fig. 1(b) (/ = 0) 
and Fig. 1(d) (/ = 12), one deduces from a continuity argument that a critical value < Ic < 
12 exists at which the two branches of the voltage nullcline intersect. 
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Figure 2: Transcritical bifurcation as the main ruler of neuronal excitability, (a) Cartoon of the 
V-nuUcline transition through a singularly perturbed transcritical bifurcation. Black circles denote stable 
fixed points, white circles unstable fixed points, (b) Continuation of the stable Ws (in green) and the 
unstable Wu (in red) manifolds of the saddle away from the singular limit {i.e. e > 0). They dictate the 
transition from the resting state (/ < Ic) to the the spiking limit cycle (7 > 7c) via a saddle-homoclinic 
bifurcation (7 — Ic). 



The critical current Ic depends on e. In the singular limit (e ~ 0) and for the corresponding 
critical current Ic — lc{0), one obtains the highly degenerate phase portrait in Figure 2a (center). 
This particular phase portrait contains a transcritical bifurcation (red circle) which is the key ruler 
of excitability. This is because, as illustrated in Figure 2b for e > 0, the convergence of solutions 
either to the resting point (/ < Ic) or to the spiking limit cycle (/ > Ic) is fully determined by the 
stable Ws and unstable Wu manifolds of the saddle point. In the singular limit shown in Figure 2a, 
these hyberbolic objects degenerate to a critical manifold that coincides with the voltage nuUcline 
near the transcritical bifurcation. It is in that sense that the X-shape of the voltage nuUcline 
completely organizes the excitability, i.e. the transition from resting state to limit cycle. 

The persistence of the manifold Wg and W^ away from the singular limit is proven by geometric 



singular perturbation (Supplementary Section S.3). The same analysis also establishes a normal 



form behavior in the neighborhood of the transcritical bifurcation: in a system of local coordinates 
centered at the bifurcation, the voltage dynamics take the simple form 



v'^ — w'^ + i + h.o.t. 



where i is a re-scaled input current and with h.o.t. referring to higher order terms in v, w, e. 

It should be emphasized that it is the same perturbation analysis that leads to the classical view 
of the Hodgkin-Huxley reduced dynamics: the transition from Figure lb left (/ = 0) to Figure lb 
left (/ — 12) involves a fold bifurcation that governs the excitability with a fold normal form 

V = v'^ — w + i + h.o.t. 

It is of interest to realize that the addition of the calcium current in the HH model unmasks 
a global view of its phase portrait that has been disregarded to date for its lack of physiological 
relevance. Figure [s] (top) shows the phase portrait of the classical reduced HH model for three 
different values of the hyperpolarizing current, revealing the transcritical singularity for the middle 
current value. The unshaded part of the first plot (and only this part of the plot) is familiar to 
most neuroscientists since the work of FitzHugh. Likewise, the conceptual sketch of the transcritical 
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3 No calcium channels (original reduced Hodgkin-Huxley model) 




Figure 3: Unfolding of the transcritical bifurcation in the global reduced Hodgkin-Huxley 
phase portrait, (a and b) Phase portraits or the original reduced HH model without (a) and with a 
calcium current (b). A constant inhibitory current of increasing amplitude (from left to right) is applied 
to the model. The transcritical bifurcation is non physiological in the classical reduced HH model (a) but 
plays an important physiological role in the presence of calcium channels (b). 

bifurcation will be familiar to all readers of basic textbooks in bifurcation analysis. For instance, 
the sketch is found in [T^ as a prototypical example of non-generic bifurcation. It is symptomatic 
that this particular example is described at length but not connected to any concrete model in 
a texbook that puts much emphasis on the relevance of bifurcation analysis in neurodynamics 
applications. As shown in Fig |3] (bottom), the missing connection is brought to life by calcium 
channels. Their particular kinetics renders the transcritical bifurcation of HH model physiological 
in the the presence of a high-conductance calcium current. 

4 Transcritical hybrid modeling of neurons 

The singular limit of planar reduced models reveals that the excitability properties of spiking 
neurons are essentially determined by a local normal form of bifurcation of the resting equilibrium. 
This property is at the core of mathematical analysis of neuronal excitability (see [H H] and the 
rich literature therein). 

In recent work, Izhikevich showed that, for computational purposes, the combination of the local 
normal form dynamics with a hybrid reset mechanism, mimicking the fast (almost discontinuous) 
spike down-stroke, is able to reproduce the behavior of a large family of neurons with a high degree 
of fidelity [131 [S] . Mimicking Izhikevich approach, we simplify the planar dynamics into the hybrid 
model: 

V = — ViP' + I if w > vt}n then (la) 

w = e{av ~ w + Wo) v -^r— c,w -^r- d (lb) 

The proposed transcritical hybrid model is highly reminiscent of the hybrid model of Izhikevich, 
but it consideraly enlarges its modeling power by including two features of importance: 

• the transcritical normal form v = — + 1 replaces the fold normal form v = — w + 1, 
in accordance with the normal form analysis of Section 3. 
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Figure 4: Schematic phase-portraits of the transcritical hybrid model for different values of / 

and wq. The v- and n-nullclines are drawn as full and dashed lines, respectively. The trajectories are drawn 
as red oriented lines. Many different phase-portraits derive from the transcritical hybrid model, including 
the one of the fold hybrid model, which only captures the shaded area. 



• the new parameter wq determines whether the intersection of the voltage and recovery null- 
clines will take place above {wq > 0) or below {wq < 0) the transcritical singularity. 

The parameter wq is a direct image of the calcium conductance: for small calcium conductances, 
the recovery variable nuUcline only intersects the upper branch of the voltage nullcline (Supplemen- 
tary Fig. ; likewise in the hybrid model when wq > 0. For high calcium conductance, the 
recovery variable nullcline intersects the lower branch of the voltage nullcline (Supplementary Fig. 
^ ; likewise in the hybrid model when wq < 0. 

Fig. 4 summarizes the four different phase portraits that derive from the transcritical hybrid 
model for different values of / and wq. For wq > 0, the model captures the classical view of the 
reduced HH model Fig. 4(bottom). For wq < 0, the model reveals the novel excitability properties 
associated to a high calcium conductance Fig. 4(top). In the Supplementary Section S.4 



we 



further investigate the hybrid transcritical phase portrait though geometrical singular perturbation 
arguments. 



5 Hybrid modeling of a thalamocortical relay neuron 

Thalamocortical (TC) relay neurons arc the input to sensory cortices. As illustrated in Figure 5, 
these neurons exhibit two distinct firing patterns: either a continuous regular spiking Fig. 5a, or 
a plateau burst spiking Fig. 5b. The switch between the two modes is regulated by prominent T- 
type calcium currents that are deinactivated by hyperpolarization, thereby modulating the resting 
membrane potential |17) . Thalamocortical neurons are widely studied in the literature, and their 
two spiking modes make them a prototypical example to illustrate the relevance of the reduced 
model. We emphasize that our objective is not a fine tuned quantitative modeling of the TC 
neuron firing pattern. Rather, we attempt to provide a qualitative picture of how the proposed 
simple hybrid dynamics permits to explain the behavior of TC neurons and, in particular, the role 
of calcium currents. 

Figure 5 compares the experimental step response of a TC neuron in vitro and the simulated 
step response of the transcritical hybrid model (Isl) , both in the low and high calcium conductance 
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Figure 5: Comparison of the experimental step response of a TC neuron in vitro [24] (top) 
and the step response of the proposed transcritical hybrid model (bottom) in low (left) and 
high calcium conductance modes (right), (a and b): Membrane potential variations of the recorded 
TC neuron over time in both conditions, (c and d) Membrane potential variations of the modeled TC 
neuron over time in both conditions. A variation of wq, which is an image of the calcium conductance, is 
sufficient to generate the switch of firing pattern physiologically observed in TC cells. 



modes. As discussed in Section [4] the small calcium conductance mode is obtained by choosing 
a positive wq, whereas the large calcium conductance mode is obtained by choosing a negative 
Wq, all the other parameters being identical in the two modes. The hybrid model reproduces the 
experimental observation: in the low-calcium mode, it responds with a slow regular train of action 
potentials; in the high-calcium mode, it responds with a long spike latency, plateau oscillations, and 
an ADP. Phase portrait analysis permits to explain the observed behavior (Supplementary Section 



S.5 1. Among other, it shows that the switch from the upper branch of the u-nullcline (low calcium 
conductance mode) to its lower branch (high calcium conductance mode) is sufficient to reproduce 
and explain the strongly different step-responses of the TC neuron. 

In order to verify the physiological consistence of the transcritical model, we further compare its 
behavior with the simulated step response of a quantitative 200-compartments model of a TC relay 
celj^ [S] in the large conductance mode (Fig. 6). For the quantitative model, we plot the trajectory 
projection on the V — d plane, where V and d denotes the somatic membrane potential and the 
activation gating variable of the somatic T-type calcium current, respectively. 

As shown on the figure, there is a striking similarity between the projection of the high dimension 
trajectory and the phase portrait of the second-order transcritical model. In both cases, the ADP 
is generated during a decrease of the activation variable, and plateau oscillations are exhibited far 
from the resting state. Moreover, the spike latency is a robust property of the transcritical model 
because the trajectory must visit the neighborhood of both the nullclines V — and w — before 
converging to the spiking limit cycle. It should be stressed that there are no comparable ways 
to reproduce this behavior in a fold hybrid model. As shown by Izhikevich [3], reproducing this 
behavior with the standard reduced HH model necessitates a non physiological alteration of the 



reset rule (see also Supplementary Sections S.6 and S.7|. This underlines the importance of the 



revisited model to capture the richness of neuron excitability. 



^Simulations were run in the Neuron environment, based on tlie configuration files freely available at 
|http: //ens ■ iaf ■ cnrs-gif ■ f r/alaln.demos .htmlj 
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Figure 6: Comparison of the step response and phase-portrait in our proposed hybrid model 
and in a 200-compartment model of TC neurons [6J, both in large calcium conductance mode. 

In the phase-portrait of the transcritical hybrid model, v- and w-nuUclines are drawn as full and dashed 
lines, respectively, and trajectories are drawn as red oriented lines. The black full line represents the v- 
nuUcline at the onset of the stimulation. The gray full line represents the w-nuUcline at the end of the burst. 
The phase-portrait of the compartmental model depicts the trajectory projection on the V — d plane, where 
V and d denotes the somatic membrane potential and the activation gating variable of the somatic T-type 
calcium current, respectively. In both cases, the ADP is generated during a decrease of the activation 
variable, and plateau oscillations are exhibited far from the resting state. 
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6 Conclusion 



The inclusion of calcium channels in Hodgkin-Huxley model has a dramatic impact on its mathe- 
matical reduction: the firing mechanisms are governed by the local normal form of a transcritical 
rather than fold bifurcation. Interestingly, it is not the phase portrait of the reduced HH model 
that is affected by calcium, but only the subregion of the plane where it is physiologically relevant. 
As a consequence, the classical FitzHugh Nagumo phase portrait is a particular (because localized) 
view of the more complete picture studied in the present paper. 

Although this enlarged phase portrait is the source of rich and diverse forms of excitability, 
its essence is captured in a simple and physiologically grounded hybrid model. The illustration 
of its modeling power on the thalamocortical neuron excitability shows the impact of revisiting 
the classical view. This illustration is just the top of the iceberg because the same principle will 
apply to many important families of neurons that are thoroughly studied and that have so far largely 
resisted reduced modeling. The proposed model will impact the understanding of excitability of e.g. 
dopaminergic, serotonergic, and subthalamic nucleus neurons, whose various firing patterns have a 
direct and critical impact in physiology and diseases, such as Parkinson's disease and depression. 
Because of its simplicity and computational efficiency, it is also an ideal candidate for physiologically 
realistic studies of high-dimensional neuronal networks. 



7 Methods 

7.1 Equation and parameters of the complete model 

The augmented HH model reads 

Hodgkin-Huxloy dynamics Calcium currents 



CV = -gKn^iV -VK)-9Nam-'h{V -VNa)-gi{V -Vl)+Iapp+Ica+Ipump 

h = a^{V){\ - n) ~ p^{V)n 
m = a.^{V){l - m) - l3m{V)m 
h = ah{V){l - h) - f3h{V)h, 

For the HH dynamics, we use the parameters of the original paper [T2]. As all other HH currents, 
the additional calcium current obeys Ohm's law 

Ica = -gCad''iV-Vca) 

where gca is the maximum calcium conductance, Vca is the calcium Nernst potential, and d is 
the calcium activation gating variable. Such kinetics fit the majority of calcium channel subtypes 
jl]. Exploiting the similarities between the potassium and calcium gating kinetics, we further 
assume that the behavior of the calcium activation gating variable d is well approximated as a 
static function of the potassium activation gating variable n. The simple choice d := n, and 
a := 3 suffices for our paper. The exponent accounts for the delayed activation of calcium currents, 
described often empirically by second to sixth powers (see [11] page 112, and references therein). 
These parameters do not reflect any precise physiological calcium current. We choose them as a 
prototypical example. The functions a^, /3a;, x = n,m, h, can be found in the paper [T^]. The value 
for the potassium Nernst potential Vk = —12 is the same as in [Hj, while the sodium Nersnt and 
the leak Nernst potential are rounded to V/va = 120 and Vi = 10.6, respectively. The values of the 
sodium gi^a = 120, potassium gx = 36, and leak gi (maximum) conductances are the same as in 
[T^ . The calcium Nernst potential is given by Vca = 150. The numerical simulations of Figure 
l(right) are obtained by picking gca — 2.7 and Ipump = —17. Parameter unchanged for Equations 
^ and Figure |3] 
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7.2 Planar reduction and phase portrait analysis 

We follow the standard reduction of the original HH model to a two dimensional system by: i) 
assuming an instantaneous sodium activation, m = rnao{V), where TOoo(^) = ctm{V) / {am(V) + 
Pm{y))', ii) exploiting the approximate linear relation, originally proposed in h + sn = c, with 
s, c ~ 1. Applying the same reduction to ^ with parameters as above, we obtain the planar system 

CV = -gKU^V - Vk) - gNamoo{Vf{0.89 - - VNa) - giiV ~ V) + hpp 

-9Can^{V - Vca) + Ipump (2a) 
fi = an{V){l ^ n) - /3n{V)n (2b) 

V^- "nuUcline" refers to the set {(V, n) G : V = 0} and similarly for other variables. 

7.3 Hybrid modeling of TC neurons 

For modeling convenience, we add a fitting parameter 6 G M in the sub-threshold (continuous) 
voltage dynamics: 

V ^ v'^ + bvw — ViP' + I 

The extra parameter does not affect the nature of the bifurcation, but tunes the slope of the 
u-nuUcline branches. In order to account for the effect of intracellular calcium variations and for 
the associated activation of calcium pump currents, we also add a slow adaptation variable z as 
follows: 

i) ^ + bvw — vo^ + I — z ifw> vtin then (3a) 

w — e{av — w + wq) v <— c, w i~ d, (3b) 

z = —e^z, z^z + dz- (3c) 

with a = 0.1, b = —3, c = 15, d = 15, e = 1, vth = 100, = 0.1, dz — 40. The low-calcium 
mode correspond to wq = 3.2, the high-calcium mode to wq = —4. The injected current / — —5 
when the stimulation is off and / = 85 when the stimulation is on. We stress two properties that 
are not captured by the reduced model ([3]): Firstly, the intracellular calcium {i.e. z) dynamics are 
generally coupled with the membrane voltage also in the subthreshold phase. Secondly, depending 
on the type of the modeled calcium current, the calcium conductance, here reflected by wq (Section 
|4]), might slowly change in response to voltage and intracellular calcium variations. Both properties 
are neglected in ([3|. 

7.4 Numerical simulations 

Numerical simulations were run with MATLAB 

(http://www.iiiathworks.com), apart from the 200-compartment model in Figures 5 and6 which 



was simulated with the NEURON software environment 

(http://www.neuron.yale.edu). The numerical bifurcation analysis of Supplementary Section S2 
has been obtained with XPP environment 



( jhttp : //www . math . pitt . edu/ ~bard/xpp/xpp . html ) 
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Supplementary material 

S.l Supplementary figures 
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Figure SI: Step responses of the HH model without (left) and with a calcium cur- 
rent (right). Time-evolution of the applied excitatory current (bottom) and of the corresponding 
membrane potential (top). 
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Figure S2: Total ionic currents for V fixed as a function of n without (left) and with 
calcium channels (right). Blank portions corresponds to the values of n where F > 0, shaded 
portions corresponds to the values of n where V <Q, and the dashed lines correspond to the values 
of n where ^ = 0. Note that the total ionic currents monotonically increase only in the absence of 
calcium channels (left). 
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Figure S3: Nullcline intersections in the reduced model (5) with calcium current for 
different values of the calcium conductance. The n-nuUcline is depicted as a dashed line, 
the y-nullchne as a solid line. The associated calcium conductance is expressed via a gray 
scale, as indicated in the figure legend. The calcium pump currents is given by Ica,pump = 
0.0, -0.74, -2.78, -15.6 for, respectively, gca = 0.0, 1.0, 1.5, 3.0. 
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S.2 Bifurcation analysis 



A bifurcation diagram of (5) with lapp as the bifurcation parameter sheds more hght on the tran- 
sition mechanism between the resting and spiking modes (Figure [l]). We use XPPAUT [1 for this 
numerical analysis. In both the /pa-on and off modes, we draw the bifurcation diagram only for 
small lapp, corresponding to the transition from resting to limit cycle oscillations (for larger lapp, 
the stable limit cycle disappears in a supercritical Andronov-Hopf bifurcation in both cases, which 
leads to a stable depolarized, i.e. high- voltage, state). 




Figure S4: One parameter bifurcation diagram of the reduced Hodgkin-Huxley model 
(5) with and without calcium current. Thin solid lines represents stable fixed points, while 
dashed lines unstable fixed points or saddle points. The thick lines labeled Vmin and Vmax represent 
the minimum and the maximum voltage of stable limit cycles, respectively, (a): Hodgkin-Huxley 
model, (b): Hodgkin-Huxley model with calcium current. Ix, with x = SNLC, AH, SH, SN , 
denotes the value of the input current for which the system undergoes the bifurcation x. More 
details are given in the text. 

Figure ^ (left) illustrates the bifurcation diagram of the original reduced Hodgkin-Huxley 
model. For low values of lapp, the unique fixed point is a stable focus that loses stability in 
a subcritical Andronov-Hopf bifurcation at lapp = I ah- Beyond the bifurcation, the trajectory 
converges to the stable spiking limit cycle. When lapp is lowered again below Isnlc, the spiking 
limit cycle disappears in a saddle-node of limit cycles, the unstable one (not drawn) emanating from 
the subcritical Andronov-Hopf bifurcation, and the trajectory relaxes back to rest. 

Figure ^ (right) illustrates the bifurcation diagram of the reduced Hodgkin-Huxley model in 
the presence of calcium channels. For lapp < Isn, a stable node (lower branch), a saddle (central 
branch), and an unstable focus (upper branch) are present, as in Figure[T}3(top right). The node and 
the saddle coalesce in a supercritical fold bifurcation at lapp = Isn, and disappear for lapp > Isn, 
letting the trajectory converge toward the stable limit cycle. The spike latency observed in the 
Jca-on configuration unmasks the ghost of this bifurcation. The stable limit cycle disappears in a 
saddle honioclinic bifurcation as lapp falls below Ish, which lets the trajectory relax back to the 
hyperpolarized state. 

The homoclinic bifurcation exhibited by the Hodgkin-Huxley model with calcium channels is 
a key mathematical difference with respect to the standard HH model. The next section unfolds 
this bifurcation by exploiting the time-scale separation between the fast voltage V and the slow 
recovery variable n. 
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S.3 Transcritical and saddle-homoclinic bifurcations of in 
the presence of calcium channels 

By exploiting the sharp time-scale separation between the membrane voltage and the recovery 
variable dynamics, in this section we assume that h — 0{e)V , where e > 0, and study the singularly 
perturbed limit e = 0. We highlight in particular the existence of a transversal V^-nuUcline self- 
intersection in the passage from Figure [IJa (top right) to Figure [T|d (bottom right) that is associated 
to a transcritical singularity (see also Fig. 2i center). We then rely on some recent results ^ to 
qualitatively study the phase portrait of (2) away from the singular limit (i.e. e > 0) and for 
different values of the input current. 



S.3.1 Transcritical bifurcation in the presence of calcium channels 

The presence of the y-nullclinc self intersection is confirmed by a closer inspection of the phase 
portrait near the nuUcline break-up, as illustrate in Figure ^ Let us discuss this transition in more 
details. The existence of the self intersection follows directly by -gp— 1 and the implicit function 
theorem [TT]. We claim that it is also transversal. Let Itc G (2.4,2.46) be the input current for 
which the ^-nuUcline has the self-intersection. Let fv{V,n, lapp) be the voltage dynamics of the 
reduced model ^ in the presence of calcium current, that is 



fviV,n,Iapp) := V\ 



=2.7, i„ 
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(S4) 



where V is defined in ([2|) and all the other parameters are unchanged. Then, we want to show that 
there exists {Vtc,ntc) G M?, such that 



fv 



dVdn 



{Vtc,ntc,Itc) 



fv{Vtc,ntc,Itc) = 
dfv, 
dV 
dfv , 
dn 

V 



dVdn 
fv 



dv? 



iVtc,ntc,Itc) = 
{Vtc,ntc,Itc) = 
{Vtc,ntcJtc) 
(Vtc^ntchc) 



< 



dVv 



{Vtc,ntc,Itc) 7^ 0, 



(S5a) 
(S5b) 

(S5c) 
(S5d) 

(S5e) 



describing a non- degenerate (transversal) self-intersection of the y-nullcline at {Vtc, ntc) (see e.g. 
[9] and Section 5.5.2, Th. 5.7 in [12]). The self-intersection satisfies the additional constraint that 
both intersecting branches are not parallel to the V axis, as implied by (55e). Condition (£ 5a I is 

do not depend on lapp, conditions (£5b) and 
< Itc varies, as in Figure ^Sla), the right and left extrema of, 



the F-nullcline equation. Noticing that and ^ 



(£5c) follows by the fact that, as /< 



app 



respectively, the left and right branches of the y-nullcline lie, by definition, on the line = 0, and, 
similarly, as lapp > Itc varies, as in Figure ^5jb), the minimum and the maximum of, respectively, 
the upper and lower branches of the l^-nuUcline lie, by definition, on the line = 0. Since at 
the intersection the four extrema coincide, conditions (£ 5b ) and (^Sc]) follow. We stress that (E 5b I 
and (£5c) define the point {Vtc, ntc) as the intersection of two lines (cf. Figure ^ 



m a unique way. 



Conditions (£5d) and (q5e| are generic and can be easily verified numerically. 

In the singular limit n = 0, the self intersection described by conditions (^5]) corresponds to 
a transcritical bifurcation (see e.g. Section Section 3.2 in [13]) of the voltage dynamics, with n 
as the bifurcation parameter. That is, as sketched in Figure [2^(center), for lapp = Itc, the two 
intersecting branches exchange their stability at n = ntc- More precisely, as n varies above and 
below the intersection point, i.e. n ^ ntc, the voltage dynamics has two fixed points, the left one 
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(a) I = 2.4 



(b) / = 2.46 



Figure S5: Phase portrait of the model for different values of the input current in the 
presence of calcium channels. The y-nuUclinc is drawn as a solid line. The dashed-lines depict 
the locus = and = 0, respectively. 



is stable, whereas the right is unstable. The stability of the fixed points follows from the fact that, 
on the left of the self-intersection, < 0, while > on the right (cf. Figure ^sj). The two 
fixed points collide in a transcritical singularity at the self-intersection. 

S.3.1.1 A normal form lemma 

In this technical section we compute a normal form of ^ associated to the self-intersection described 
by conditions (^. 
Give e > 0, let 

eg{V,n) := n, (S6) 

where n is defined in ([2]), and 

90 g{Vtc,ntc), (S7) 



where Vtc and rite satisfy the defining conditions (q5b|-(£5b|. The form is just a rescaling 
(through e) of the recovery variable dynamics that highlights the timescale separation between V 
and n. The following lemma is an application of Lemma 2.1 in |S] to the reduced HH model ^ 
with calcium current. 

Lemma 1 Suppose that gg < 0. Then there exists an ajjine change of coordinates and a rescaling 
of e and I that transforms the dynamical system 

V = fv{V, ) (S8a) 

n = eg{V,n), (SSb) 

where fy and g are respectively defined in ( ^J) and ( into the dynamical system 

V = v'^ - w'^ + I + hi{v,w,e) (S9a) 
w = e(-l + /i2(f,w,e)), (S9b) 



where hi{v, w, e) ~ 0{v^, v'^w, vw^ , , ev, ew, e^) and h2{v, w, e) = 0{v, w, e). 
Pro 

7^ 



Proof Let a (Fte,ntc,/tc), /? ^§^{Vtc,ntc,Itc), 1 ■= (Vtc, rite, he), and A := 

. From (5) and PI Lemma 2.1], it follows that, for / = Itc, the affine change of variable 
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(c) I = IsH~ 2.497 (d) / = 2.52 > Ish 



Figure S6: Stable and unstable manifolds of the saddle point as lapp varies. (a,b) For 
lapp < IsH, the left branch of the unstable manifold of the saddle Wu (in violet) returns (after a 
spike generation) in the vicinity of the saddle point on the left of the stable manifold (in green) 
and toward rest, (c) At lapp = Ish, there is a homoclinic orbit F (in red), connecting the stable 
and the unstable manifolds, (d) For lapp > Ish, the left branch of the unstable manifold returns 
in the vicinity of the saddle on the right of the stable manifold and converges toward the spiking 
limit cycle attractor. 

V = a{V — Vtc) + I3{n — ntc), w — \/ ii"^ ~ ^ain — rite), transforms after a suitable rescaling of 
e, into the equation 

w = — + Ae + hi{v, w, e) 
ih = e(— 1 + h2{v,w,e)). 

Noticing that /y is afhne in the input current, the extra term a{I — Itc) must be added to ii in the 
case / 7^ Itc- The result follows by defining the rescaled input current J := Ae + a{I — Itc)- D 

Note that for the dynamical system (^9]) with / = 0, the self intersection of the y-nuUcline 
discussed in Section fS . 3 . 1 1 becomes evident, with the two intersecting branches given hy w = ±v. 

S.3.2 Singularly perturbed saddle-homoclinic bifurcation in the presence 
of calcium channels 

Figure ^ provides a close look of the saddle-homoclinic bifurcation in the reduced HH model ^ 
with calcium current. In order to rigorously prove the existence of this bifurcation, we rely on 
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the normal form introduced in Lemma [T] and exploit the timescale separation between v and w 
through geometrical singular perturbations theory (GSPT). To this aim we must briefly recall some 
basic of GSPT, relying on (^9]) as an explicit example. The interested reader will find in [7] an 
excellent introduction to the topic, and in [101 [HI [S] some recent extensions on which we rely for the 
forthcoming analysis. 

The time rescaling r :— et transforms (q9|) into the equivalent system 



ev 



hi{v, w, e) 



w = {-I + h2{v,w,e)), 



(SlOa) 
(SlOb) 



which describes the dynamics (q9^ 
to as the singular limit, one obtains from 
dynamics 

= 1 

w 



in the slow timescale r. In the limit e = 0, commonly referred 
91) and (£ 10 ) two new dynamical systems: the reduced 



v'^ — w'^ + I + hi {v, w, e) 
(-1 + h2{v,w,e)), 



evolves in the slow timescale r, while the layer dynamics 

V ~ v'^ — w'^ + I + hi{v,w, e) 
w = 0, 



(Slla) 
(Sllb) 

(S12a) 
(S12b) 



evolves in the fast timescale t. Figure q7|[a) depicts the fast-slow dynamics associated to (£ 11 )- 
(£12). The main idea behind GSPT is to combine the analysis of the reduced and layer dynamics 
to derive conclusion on the behavior of the nominal system, i.e. with e > 0. 



The reduced dynamics (£11 1 is a dynamical system on the set 5*0 := {{v,w) G 



+ / + hi{v, w, e) = 0}, usua lly c alled the critical manifold. The points in So are indeed critical 



points of the layer dynamics (£12). More precisely, portions of 5*0 on which is non- vanishing 



are normally hyperbolic invariant manifolds of equilibria of the layer dynamics, whose stability 



is determined by the sign of 



dv 

dV 



Conversely, points in Sq where 



dV 

av 



constitute degenerate 



equilibria. In particular, the layer dynamics (£12) exhibits, for / = 0, two degenerate equilibria. As 
depicted in Figure ^7]^ a), they are given by the self-intersection of the V^-nuUcline, corresponding to 
a transcritical singularity (Section S.3.1), and by the fold singularity at the maximum of the upper 
branch of the V^-nuUcline. 

The basic result of GSPT, due to Fennichel is that, for e sufficiently small, non-degenerate 
portions of Sq persist as nearby normally hyperbolic locally invariant manifolds S'^ of (^8]). More 
precisely, the slow manifold S'^ lies in a neighborhood of 5*0 of radius 0(e). The dynamics on S'' 



is a small perturbation of the reduced dynamics (£11). We point out that S"^ may not be unique, 
but is determined only up to 0(e~'^/'^), for some c > 0. That is, two different choices of S'^ are 
exponentially close (in e) one to the other. Since the presented results are independent of the 
particular S"^ considered, we let this choice be arbitrary. The trajectories of the layer dynamics 
perturb to a stable and an unstable invariant foliations with basis S'^. 

The analysis near degenerate points is more delicate. Only recently some works have treated 
this problem in its full generality for different types of singularities (TUIOIH]. Figure ^(b),(c),(d) 
sketch the extension of the attractive slow manifold S^ after the fold point, and the three possible 
ways in which S^ and the repelling slow manifold S^ can continue after the trascritical singularity, 
depending on the injected current. 

The result depicted in Figure ^relies on the following theorem, adapted from [5]. 

Let A := {{v,w) € : w_ < w < w+, w — p}, be the section depicted in Figure ^ where 
p < and IpI is sufficiently small, and V-,v+ are such that A Ci S~ 7^ 0. For a given e > 0, 
let qa,e ■= A n SI and Qr.t :— A Ci S^ be the intersections, whenever they exist, of respectively 
the attractive and repelling invariant submanifolds S^ and S^ with the section A. The following 
theorem reformulates in a compact way the discussion contained in Remark 2.2 and Section 3 of 
Pj^for systems with inputs of the form (^9|. 



The first author is thankful to Prof. Szmolyan for his useful comments. 
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Figure S7: Phase-portrait of (^9]) for different values of e and /. (a): fast-slow dynamics of 
(^9| for e ~ I ~ 0. The attractive (resp. repelling) branches of the critical manifold 5*0 above/below 
the transcritical singularity are denoted by /S~ (resp. S:^/S~). (b): Continuation of the slow 
attractive and repelling manifolds for e > and / < Ic{y/e), where Ic{\/e-) is defined as in 
Theorem [l] (c): Continuation of and for e > and / — Ic{\/f)- (d): Continuation of 5*^ and 
SI for e > and / > /c(%A). 
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Theorem 1 (Adapted from [9]) Consider the system (^^. Then there exists eq > and a 
smooth function Ic{\/e), defined on [0, Eq] and satisfying /c(0) = 0, such that, for all e G (0, Eq], the 
following assertions hold 

1. qa,e = qr.c if and only if I ^ IciV^) 

2. there exists an open interval A 3 Ic{\/e), such that, for all I ^ A, it holds that A H 5^ 7^ 0, 
A n S"^ 7^ 0, and 

d , ^ „ 

Figure ^illustrates this result. 

Remark 1 The function IdV^) is related to the function Ac(-\/e) defined in [51 Remark 2.2] by 
'■— eAc(v^). Similarly, given e > 0, the parameter / appearing in Theorem [l] is just the 
re-scaling / = eA of the parameter A appearing in [9l Remark 2.2 and Sections 3]. 

Theorem [T] implies the existence of the saddle-homoclinic bifurcation in the reduced Hodgkin- 
Huxley model ^ with calcium current. As stressed by Figure ^b,c,d), the slow attractive 
(resp. repelling S^) manifold coincides with the unstable Wu (resp. stable Ws) manifold of the 
saddle point, as it can be proved via qualitative arguments. Thus, for / < /c(-/e), the unstable 
manifold Wu continues after the transcritical singularity on the left of Wg, toward the stable node. 
See Figure ^a,b) and Figure ^b). For I = /c(-\/e), Wu extends after the transcritical point to 
Wg, forming the saddle- homocHnic trajectory, as depicted in Figure ^c) and Figure ^c). For 
I > Ici^/e), the unstable manifold of the saddle Wu continues after the transcritical singularity 
on the right of Wg , and spirals toward an exponentially stable limit cycle, whose existence can be 
proved with similar GSPT arguments (see for instance [lO]). This situation is the one depicted in 
Figure ^d) and Figure ^d). 



S.4 Hybrid singularly perturbed saddle-homoclinic bifurca- 
tion 

The saddle-homoclinic bifurcation analysis provided in Section [S. 3. 2| for the reduced calcium-gated 
HH model naturally extends to the hybrid dynamics ([T]) with wq < 0. Indeed, by construction, if 
Wo < 0, this model can be transformed in the normal form (^9| derived in Lemma [ij and Theorem 
[ijappHes directly. 

Similarly to the derivation in Section S.3.2 we can associate to ([!]) two new dynamical systems, 



describing its singular dynamics in the fast timescale t and in the slow timescale r := et, respectively. 
More precisely, the hybrid reduced dynamics 

= ~ + I if V > Vfh, then 

w = av — w + wq V c, u ^ d, 

evolves in the slow time scale r, while the hybrid layer dynamics 

V = v'^ ~ w"^ + I if w > vth, then 

w = V -i^ c, u ^ d, 

evolves in the fast time-scale t. Figure ^Sja) depicts the associated slow-fast hybrid dynamics for 
/ = 0, a e (0, 1), and wq < 0. 

The analysis of the non-singular limit follows the same line as the analysis developed in Section 
|S.3.2| for the continuous time case. The only difference is that the return mechanism, provided in 
the continuous time case by the right attractive branch of the critical manifold Sq, is now replaced 
by the hybrid reset. The result is summarized in Figure ^(b),(c),(d). The slow attractive manifold 

is chosen as the continuation of the trajectory starting at the reset point. In this way, it also 
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Figure S8: Phase-portrait of the hybrid dynamics ([T]) for different values of e and /. 

(a): Fast-slow dynamics for e = / = 0. The attractive (resp. repelling) branches of the critical 
manifold 5*0 above/below the transcritical singularity are denoted by S^/S~ (resp. S^/S^). (b): 
Continuation of the slow manifolds for e > and / < Idy/e), where /c(-\/e) is defined as in Theorem 
[l] (c): Continuation of the slow manifolds for e > and / = Ic{\f(-)- (d): Continuation of the slow 
manifolds for e > and / > Ic{\/e). 
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coincides with the image of the unstable manifold of the saddle Wu through the hybrid reset map. 
As in the continuous time case, the stable manifold of the saddle Wg can be shown to coincide with 
the slow repelling manifold 5*^. Let IdV^) be defined as in Theorem[l] For / < Ic{^/e), the unstable 
manifold W„ is brought back on the left of Wg , and the resting state is globally stable, as in Figure 
^b). At / = Ic{y/e), the hybrid reset connects Wu and Wg, corresponding to a hybrid homoclinic 
bifurcation. The associated phase-portrait is shown in Figure ^Sj^c). Finally, for / > /c(\/e), the 
unstable manifold Wu is brought by the reset mechanism on the right of Wg and directly into the 
newborn hybrid limit cycle attractor, as in Figure ^8j[d). 

We emphasize that the existence of the hybrid homoclinic bifurcation and the associated critical 
value Ic are independent of the reset point (c, d), provided that the trajectory starting from (c, d) 
be attracted toward the left upper branch of the critical manifold, as in Figure ^ As discussed 
in Section |S.3.2[ under this condition, two different choices of the reset point are associated to two 
slow manifolds that are near, for some c > 0. By the result in Theorem [l] the two 

values of the critical current Ic for which the slow attractive manifold extends to the slow repelling 
manifold are again near. This also ensures that the value for which the hybrid homoclinic 

bifurcation happens is independent of the reset point, modulo variations that are 0(e~'^/'^). The 
same robustness properties are not shared by fold hybrid models (see Section S.7 below). 



S.5 Phase-plane analysis of the transcritical hybrid TC neu- 
ron model 

Fig. ^ shows the time-course of the qualitative transcritical hybrid modeling of thalamocortical 
relay (TC) neurons ([S]) in the low and high calcium conductance modes, as well as the corresponding 
phase portraits. The analysis of the phase-portraits gives insights on the mechanisms by which TC 
cells exhibit tonic firing or bursting when submitted to a similar step input, according to the initial 
resting potential. 

When wq > (Fig. ^ left), which corresponds to a small calcium conductance (T-type calcium 
channels are inactivated), the hyperpolarized state belongs to the upper branch of the w-nullcline. 
Application of a depolarizing current step lifts the voltage nuUcline above the resting state, thus 
generating a transient non-delayed action potential (marked with a * in Fig. ^9jD, left). Note 
that, when the hyperpolarized state belongs to the upper branch of the w-nuUcline, no plateau 
oscillations are possible (Fig. ^9|d, left). Furthermore, the relaxation toward the hyperpolarized 
state is necessarily monotone (i.e. no ADP), as stressed in Fig. ^9j:, left. 

At the generation of the first spike, a certain amount of calcium ions enter the cell due to the 
presence of high voltage activated calcium currents. The subsequent activation of outward calcium 
pump currentfj^transiently hyperpolarizes the cell. As the intracellular calcium is expelled (Fig. [9]:, 
left), calcium pump currents slowly deactivate and the cell slowly depolarizes (interspike period), 
until the spiking threshold is reached and a new action potential is fired. 

When Wo < (Fig. ^ right), which corresponds to a large calcium conductance (T-type 
calcium channels are available), the hyperpolarized state belongs to the lower branch of the v- 
nuUcline, and is more hyperpolarized than for wq > 0, as in experiments. When a depolarizing 
current step is applied, this branch falls below the w-nullcline (Fig. ^9]b, right). In order to 
generate the first spike, the state travels in the narrow region between the two nuUclines, resulting in 
a pronounced latency. Furthermore, the relative position of the hyperpolarized state (lower branch) 
with respect to (hybrid) spiking cycle (upper branch) clearly explains the generation mechanism of 
plateau oscillations. These high frequency plateau oscillations (burst) continue until the intracellular 
calcium, which enters at each action potential, and the associated pumps current {i.e. z) are 
sufficiently large. Plateau oscillations then terminate in a (hybrid) saddle-homoclinic bifurcation 
(Fig. ^9]:, right). At the end of the burst, the system converges toward the hyperpolarized state 
following the left branch of the w-nuUcline, thus generating a marked ADP at the passage near the 
nullcline's funnel. The subsequent slow phase is mainly ruled by the variations of the intracellular 

''This activation is modeled by the hybrid reset z <— z + d^. 
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calcium. With the adopted simple dynamics it consists in a slow depolarization that follows the 
decrease of the intracellular calcium (Fig. right). A finer and more physiological modeling 

of the intracellular calcium dynamics could reproduce in vitro recordings with a higher degree of 
fidelity. 



S.6 Modeling TC neurons with a fold hybrid model 

This section reproduces the modeling of TC neuron with the fold hybrid model. Such a model 
was recently proposed by Izhikevich and used in a large-scale model of the mammalian thalamic 
system [B]. The model succeeds in reproducing the firing pattern of TC neurons by modifying 



the reset mechanism. Fig. £ 10 shows that the voltage time-course reproduces with fidelity the 
three hallmarks of large calcium conductances, that is (A) spike latency, (B) plateau oscillations, 
and (C) ADP. However, the phase portrait illustrates that the modified u-reset law has no obvious 
physiological interpretation. As a consequence, plateau oscillations are not endogenously generated 
in this model. Furthermore, as opposed to both the quantitative model and the proposed hybrid 
model, the ADP is generated following an increase of the recovery variable, when the reset point 
crosses the stable manifold of the saddle. 



S.7 Robust ADP generation 

For a neuron model to be biologically relevant, it should be robust to exogenous disturbances (small 
synaptic inputs, thermal noise, etc.). The firing pattern, in particular, should remain unchanged. 
Figure compares the perturbation robustness of three TC neuron models to small current 
impulses. It suggests that the fold hybrid model is less robust than the transcritical model, because 
a tiny pulse is sufficient to generate an extra action potential at the ADP apex. 

The difference in robustness is explained by the different ADP generation mechanisms, as illus- 



trated in Figure £ 12 In the fold model, ADPs are generated when trajectories cross the w-nuUcline 



from below (cf. Figures £ 10 and £ 12). The absence of any robust attractor in the ADP generation 
region makes the ADP height and shape heavily dependent on the exact reset point. Moreover, 
when small current pulses are applied, the ADP generation is disrupted, and the model fires an 
extra (non-physiological) spike. 

Conversely, ADP generation in the transcritical hybrid model is robustly governed by the attrac- 
tor SI that steers the trajectories through the ADP apex and toward the resting point. That is the 
reason why the ADP height and shape barely depend on chosen reset point. At the same time, the 
persistence to small perturbations of this invariant manifold [5] ensures, as required in biologically 
meaningful conditions, the robustness of the ADP generation mechanism to small inputs. 
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Figure S9: Transcritical hybrid modeling of TC cell ([s]) in the low {wq — 3.2, left) 
and high calcium conductance modes (wq = 4, right), (a) Time-courses of the model in 
both conditions, (b and c) Phase-portrait of the transcritical hybrid model in both conditions. 
Trajectories are depicted as solid oriented red lines. The reset point is depicted as a square □, 
while the (instantaneous) hyperpolarized state as a filled circle •. The w-nullcline is depicted as a 
dashed line. In (b), the gray (black) thin solid line is the u-nuUcline when the current step in off 
(on). In (c), the w-nuUcline is depicted as gray thin lines of different darkness. As sketched in the 
figure, light gray correspond to large calcium concentrations, whereas dark gray to small. 
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Figure SIO: Phase portrait and voltage time course (inset) of the Izhikevich model of 
TC neurons. The trajectory is depicted as a solid oriented red line. The w-nulcline is depicted 
a solid line, the w-nuUcline as a dashed line. The point line sketches the w-rest law. The stable 
manifold of the saddle point is depicted as a green line. Parameters as in [4, Sections 8.3.1]. 
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Figure Sll: Nominal step response (left) and step response in the presence of small 
current pulses in the 200 compartments TC neuron model (a), the proposed hybrid 
model (b), and the Izhikevich model of TC neuron ([5|, Section 8.3.1]) (c). 




Figure S12: Comparison of the ADP generation mechanisms in the fold (left) and in the 
transcritical hybrid models (right). The stable manifold of the saddle (x) is depicted in green. 
In the Izhikevich model ADPs are generated by sliding near the stable manifold of the saddle and 
crossing the w-nuUcline from below. In our proposed hybrid model, ADP are robustly generated 
along the attractor S^. 
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